interactive HALO tracks¶
This notebook shows how to plot the flight tracks of the HALO aircraft. First, we’ll have only a quick look and later on we do some more work in order to get the data onto an interactive map.
Preparations¶
First, we’ll import pylab and the EUREC4A data catalog.
%pylab inline
import eurec4a
cat = eurec4a.get_intake_catalog()
Populating the interactive namespace from numpy and matplotlib
Restructuring the dataset¶
When inspecting the contents of the currently available BAHAMAS quicklook dataset, we see that it does not handle coordinates according to the ideas of CF conventions, which makes using the data a lot more difficult.
ds = cat.HALO.BAHAMAS.QL['HALO-0122'].to_dask()
ds
<xarray.Dataset>
Dimensions: (tid: 331760)
Dimensions without coordinates: tid
Data variables: (12/46)
TIME (tid) datetime64[ns] ...
IGI_RMSX (tid) float32 ...
IGI_RMSY (tid) float32 ...
IGI_RMSZ (tid) float32 ...
IRS_ALT (tid) float32 ...
IRS_ATA (tid) float32 ...
... ...
TS (tid) float32 ...
RELHUM (tid) float32 ...
SOURCE (tid) float32 ...
WS (tid) float32 ...
MIXRATIOV (tid) float32 ...
MIXRATIO (tid) float32 ...
Attributes: (12/20)
title: Quality controlled data of DLR measurment flight...
aircraft: D-ADLR
project: EUREC4A
mission: EUREC4A
ProjectInvestigator: Stevens
flightname: adlr_20200122a
... ...
TimeInterval: 14:57:36 - 24:10:32
comment: Humidity is measured with SHARC and/or HUMICAP ...
platform: HALO
instrument: BAHAMAS
Variablelist: TIME IGI_RMSX IGI_RMSY IGI_RMSZ IRS_ALT IRS_ATA...
EXTRA_DIMENSION.freqid: 10- tid: 331760
- TIME(tid)datetime64[ns]...
- long_name :
- time in seconds since 1970-1-1 (UTC)
- standard_name :
- strptime_format :
- double
[331760 values with dtype=datetime64[ns]]
- IGI_RMSX(tid)float32...
- units :
- m
- long_name :
- RMS/Position North
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IGI_RMSY(tid)float32...
- units :
- m
- long_name :
- RMS/Position East
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IGI_RMSZ(tid)float32...
- units :
- m
- long_name :
- RMS/Position Up
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_ALT(tid)float32...
- units :
- m
- long_name :
- WGS84 Datum/Elliptical Height
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_ATA(tid)float32...
- units :
- deg+
- long_name :
- IGI actual track angle
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_AXB(tid)float32...
- units :
- m/s^2
- long_name :
- IGI body x-axis acceleration
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_AYB(tid)float32...
- units :
- m/s^2
- long_name :
- IGI body y-axis acceleration
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_AZB(tid)float32...
- units :
- m/s^2
- long_name :
- IGI body z-axis acceleration
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_AZG(tid)float32...
- units :
- m/s^2
- long_name :
- IGI vertical acceleration
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_EWV(tid)float32...
- units :
- m/s
- long_name :
- velocity/East
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_GS(tid)float32...
- units :
- m/s
- long_name :
- IRS Groundspeed from corrected IGI data
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_HDG(tid)float32...
- units :
- deg+
- long_name :
- true heading
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_LAT(tid)float64...
- units :
- deg
- long_name :
- WGS84 Datum/Latitude
- standard_name :
- strptime_format :
- double
[331760 values with dtype=float64]
- IRS_LON(tid)float64...
- units :
- deg
- long_name :
- WGS84 Datum/Longitude
- standard_name :
- strptime_format :
- double
[331760 values with dtype=float64]
- IRS_NSV(tid)float32...
- units :
- m/s
- long_name :
- velocity/North
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_P(tid)float32...
- units :
- deg/s
- long_name :
- IGI body roll rate
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_PHI(tid)float32...
- units :
- deg
- long_name :
- roll angle
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_Q(tid)float32...
- units :
- deg/s
- long_name :
- IGI body pitch rate
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_R(tid)float32...
- units :
- deg/s
- long_name :
- IGI body yaw rate
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_THE(tid)float32...
- units :
- deg
- long_name :
- pitch angle
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- IRS_VV(tid)float32...
- units :
- m/s
- long_name :
- velocity/Up
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- ABSHUM(tid)float32...
- units :
- g/m^3
- long_name :
- Absolute Humidity
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- ALPHA(tid)float32...
- units :
- deg
- long_name :
- Angle of Attack
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- BETA(tid)float32...
- units :
- deg
- long_name :
- Angle of Sideslip
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- H(tid)float32...
- units :
- m
- long_name :
- Height above Sea Level calculated from met. data (fast)
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- HP(tid)float32...
- units :
- m
- long_name :
- Pressure Altitude
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- MC(tid)float32...
- units :
- Ma
- long_name :
- Machnumber
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- PS(tid)float32...
- units :
- hPa
- long_name :
- Static Pressure (from NB_PSIA)
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- QC(tid)float32...
- units :
- hPa
- long_name :
- Dynamic Pressure
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- TAS(tid)float32...
- units :
- m/s
- long_name :
- Calculated True Airspeed
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- TAT(tid)float32...
- units :
- K
- long_name :
- Total Air Temperature (from BDY-TTQ, anti-ice-corrected)
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- TD(tid)float32...
- units :
- K
- long_name :
- Dewpoint Temperature
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- THETA(tid)float32...
- units :
- K
- long_name :
- Potential Temperature
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- THETA_V(tid)float32...
- units :
- K
- long_name :
- Virtual Potential Temperature
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- TV(tid)float32...
- units :
- K
- long_name :
- Virtual Temperature
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- U(tid)float32...
- units :
- m/s
- long_name :
- Wind Vector East Component
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- V(tid)float32...
- units :
- m/s
- long_name :
- Wind Vector North Component
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- W(tid)float32...
- units :
- m/s
- long_name :
- Wind Vector Vertical Component
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- WA(tid)float32...
- units :
- deg+
- long_name :
- Horizontal Wind Direction
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- TS(tid)float32...
- units :
- K
- long_name :
- Static Air Temperature (from BDY-TTQ, anti-ice-corrected)
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- RELHUM(tid)float32...
- units :
- %
- long_name :
- Relative Humidity with Respect to Water
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- SOURCE(tid)float32...
- units :
- long_name :
- Humidity Data Source (1:Vaisala/2:SHARC_longpath/3:SHARC_shortpath)
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- WS(tid)float32...
- units :
- m/s
- long_name :
- Horizontal Windspeed
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- MIXRATIOV(tid)float32...
- units :
- micromol/mol
- long_name :
- H2O Volume Mixing Ratio
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- MIXRATIO(tid)float32...
- units :
- g/kg
- long_name :
- H2O Mass Mixing Ratio
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
- title :
- Quality controlled data of DLR measurment flight adlr_20200122a
- aircraft :
- D-ADLR
- project :
- EUREC4A
- mission :
- EUREC4A
- ProjectInvestigator :
- Stevens
- flightname :
- adlr_20200122a
- flightnumber :
- 471
- source :
- Flight Facility DLR Oberpfaffenhofen
- system :
- Data from D-ADLR measurement system
- institute :
- Flight Facility DLR Oberpfaffenhofen
- contact :
- A.Giez, V.Dreiling, M.Zoeger, Ch. Mallaun; email: andreas.giez@dlr.de
- date :
- 2020-1-22
- date_created :
- 2020-1-22
- date_last_revised :
- 2020-1-25
- TimeInterval :
- 14:57:36 - 24:10:32
- comment :
- Humidity is measured with SHARC and/or HUMICAP VAISALA instruments with emphasis on the SHARC TDL instrument. The humidity source is defined in the variable \"source\": 0: no data 1: VAISALA 2: SHARC short path 3: SHARC long path. Within clouds, the humidity measurement can be affected by evaporation of cloud particles resulting in rel. hum exceeding 100%. Sensor wetting within clouds can lead to wrong temperature measurements. This error can be up to a few Kelvin. The misalignment of IGI platform compared to acft-IRS is found as (rad): pitch= 0.00834372, roll= -0.00012623
- platform :
- HALO
- instrument :
- BAHAMAS
- Variablelist :
- TIME IGI_RMSX IGI_RMSY IGI_RMSZ IRS_ALT IRS_ATA IRS_AXB IRS_AYB IRS_AZB IRS_AZG IRS_EWV IRS_GS IRS_HDG IRS_LAT IRS_LON IRS_NSV IRS_P IRS_PHI IRS_Q IRS_R IRS_THE IRS_VV ABSHUM ALPHA BETA H HP MC PS QC TAS TAT TD THETA THETA_V TV U V W WA TS RELHUM SOURCE WS MIXRATIOV MIXRATIO
- EXTRA_DIMENSION.freqid :
- 10
In order to fix that a bit, we could come up with a little function which transforms the dataset we get into a dataset we would like to have. Of course the proper way of handling this would be to actually fix the dataset, but this way of doing it may still be illustrative. Another advantage of this preparatory step is that amount of data in the dataset can be reduced, such that when we later load the datasets early, less data will be fetched.
def fix_halo_ql(ds):
import xarray as xr
ds = ds.rename({"tid":"time"})
return xr.Dataset({
"time": ds.TIME,
"lat": ds.IRS_LAT,
"lon": ds.IRS_LON,
"alt": ds.IRS_ALT,
})
This at least has a time coordinate.
dsfixed = fix_halo_ql(ds)
dsfixed
<xarray.Dataset>
Dimensions: (time: 331760)
Coordinates:
* time (time) datetime64[ns] 2020-01-22T14:57:36 ... 2020-01-23T00:10:3...
Data variables:
lat (time) float64 ...
lon (time) float64 ...
alt (time) float32 ...- time: 331760
- time(time)datetime64[ns]2020-01-22T14:57:36 ... 2020-01-...
- long_name :
- time in seconds since 1970-1-1 (UTC)
- standard_name :
- strptime_format :
- double
array(['2020-01-22T14:57:36.000000000', '2020-01-22T14:57:36.100000000', '2020-01-22T14:57:36.200000000', ..., '2020-01-23T00:10:31.700000000', '2020-01-23T00:10:31.800000000', '2020-01-23T00:10:31.900000000'], dtype='datetime64[ns]')
- lat(time)float64...
- units :
- deg
- long_name :
- WGS84 Datum/Latitude
- standard_name :
- strptime_format :
- double
[331760 values with dtype=float64]
- lon(time)float64...
- units :
- deg
- long_name :
- WGS84 Datum/Longitude
- standard_name :
- strptime_format :
- double
[331760 values with dtype=float64]
- alt(time)float32...
- units :
- m
- long_name :
- WGS84 Datum/Elliptical Height
- standard_name :
- strptime_format :
- float
[331760 values with dtype=float32]
Just to have a better visual impression, we can create a quick overview plot of the data:
plt.plot(dsfixed.lon, dsfixed.lat)
plt.show()
Reducing the size of the dataset¶
Later on, we want to plot all flights on an interactive map. Currently the dataset is rather large as the aircraft location has been recorded continuously at a high data rate. While this is good for quantitative analysis, this leads to poor interactive performance. So before going further, it is a good idea to reduce the amount of data while keeping the visual impression.
A possible idea to reduce the amount of required data is that plotting a line already does linear interpolation between two coordinate points. So if we would remove those points which are close to the linear interpolation between its neighboring points, the visual impression will stay almost the same. This idea has already been stated by Ramer Douglas and Peucker and is illustrated at Wikipedia. While the algorithm is not hard to write, it is difficult to do it efficiently in python. Thus, I’ve skipped it and use the simplification library in stead.
Note
Many algorithms for shape processing are contained in the beautiful shapely library. In general, I’d recommend using that library, but it requies the GEOS library which is a bit tricky to install.
from simplification.cutil import simplify_coords_idx
def simplify_dataset(ds, tolerance):
indices_to_take = simplify_coords_idx(np.stack([ds.lat.values, ds.lon.values], axis=1), tolerance)
return ds.isel(time=indices_to_take)
We can now use that algorithm to generate a simplified version of the dataset with a tolerance of \(10^{-5}\) degrees, which otherwise looks the same to the previous version.
dssimplified = simplify_dataset(dsfixed, 1e-5)
dssimplified
<xarray.Dataset>
Dimensions: (time: 8225)
Coordinates:
* time (time) datetime64[ns] 2020-01-22T14:57:36 ... 2020-01-23T00:10:3...
Data variables:
lat (time) float64 13.08 13.08 13.08 13.08 ... 13.07 13.07 13.07 13.07
lon (time) float64 -59.49 -59.48 -59.48 -59.47 ... -59.52 -59.51 -59.5
alt (time) float32 7.656 60.04 130.9 194.4 ... 139.0 79.67 34.11 5.528- time: 8225
- time(time)datetime64[ns]2020-01-22T14:57:36 ... 2020-01-...
- long_name :
- time in seconds since 1970-1-1 (UTC)
- standard_name :
- strptime_format :
- double
array(['2020-01-22T14:57:36.000000000', '2020-01-22T14:57:44.000000000', '2020-01-22T14:57:49.600000000', ..., '2020-01-23T00:10:01.700000000', '2020-01-23T00:10:16.700000000', '2020-01-23T00:10:31.900000000'], dtype='datetime64[ns]')
- lat(time)float6413.08 13.08 13.08 ... 13.07 13.07
- units :
- deg
- long_name :
- WGS84 Datum/Latitude
- standard_name :
- strptime_format :
- double
array([13.075457, 13.077002, 13.078168, ..., 13.068407, 13.070375, 13.072349])
- lon(time)float64-59.49 -59.48 ... -59.51 -59.5
- units :
- deg
- long_name :
- WGS84 Datum/Longitude
- standard_name :
- strptime_format :
- double
array([-59.488775, -59.482688, -59.478265, ..., -59.517247, -59.509427, -59.501452]) - alt(time)float32...
- units :
- m
- long_name :
- WGS84 Datum/Elliptical Height
- standard_name :
- strptime_format :
- float
array([ 7.655899, 60.038113, 130.8742 , ..., 79.667366, 34.112293, 5.527757], dtype=float32)
We can now compare those two tracks side by side while keeping a look at the number of data points required.
fig, (ax1, ax2) = subplots(1, 2, figsize=(15,4))
ax1.plot(dsfixed.lon, dsfixed.lat)
ax1.set_title(f"{len(dsfixed.time)} data points")
ax2.plot(dssimplified.lon, dssimplified.lat)
ax2.set_title(f"{len(dssimplified.time)} data points")
plt.show()
ratio = len(dssimplified.time) / len(dsfixed.time)
print(f"compression ratio: {ratio*100:.2f} %")
compression ratio: 2.48 %
The dataset size has been substantially reduced while the visual impression stayed the same.
First interactive map¶
In order to show the map a little bit more interactively, we use the ipyleaflet library which creates a bridge between ipython and the Leaflet JavaScript library.
import ipyleaflet
As we will need to convert many tracks to ipyleaflet layers later on, the easiest is to create a little function for that purpose right away:
def track2layer(track, color="green", name=""):
return ipyleaflet.Polyline(
locations=np.stack([track.lat.values, track.lon.values], axis=1).tolist(),
color=color,
fill=False,
weight=2,
name=name
)
With the help of that little function, creating a map is now like a breeze:
testmap = ipyleaflet.Map(center=(13.3, -57), zoom=7)
testmap.add_layer(track2layer(dssimplified))
display(testmap)
All in one¶
Let’s see if we can add all the flights and provide a layer switcher so that we can have a look at all flights individually. We’ll start by loading and simplifying all out track data into a local dictionary. Requesting the datasets from the server is IO bound and thus can simply be accellerated using a ThreadPool:
from multiprocessing.pool import ThreadPool
pool = ThreadPool(20)
def get_dataset(flight_id):
ds = cat.HALO.BAHAMAS.QL[flight_id].to_dask()
return flight_id, fix_halo_ql(ds).load()
full_tracks = dict(pool.map(get_dataset, cat.HALO.BAHAMAS.QL))
We still have to simplify the dataset, which is done here:
tracks = {flight_id: simplify_dataset(ds, 1e-5)
for flight_id, ds in full_tracks.items()}
Let’s also quickly grab some colors from a matplotlib colorbar, such that we can show all tracks in individual colors:
colors = [matplotlib.colors.to_hex(c)
for c in plt.cm.inferno(np.linspace(0, 1, len(tracks)))]
We can now start with a new empty map. Let’s also have a different basemap.
m = ipyleaflet.Map(
basemap=ipyleaflet.basemaps.Esri.NatGeoWorldMap,
center=(13.3, -57), zoom=7
)
We’ll add all the tracks as individual layers to the map
for (flight_id, track), color in zip(tracks.items(), colors):
m.add_layer(track2layer(track, color, flight_id))
and add a scale, a legend, layer controls and a full screen button to the map and show it. If you want to zoom in, you can for example shift-click and drag a rectangle over the area you want to zoom in more closely.
m.add_control(ipyleaflet.ScaleControl(position='bottomleft'))
m.add_control(ipyleaflet.LegendControl(dict(zip(tracks, colors)),
name="Flights",
position="bottomright"))
m.add_control(ipyleaflet.LayersControl(position='topright'))
m.add_control(ipyleaflet.FullScreenControl())
display(m)